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Abstract 

We study the flow of an incompressible Uquid film down a wavy incline. Applying a Galerkin 
method with only one ansatz function to the Navier-Stokes equations we derive a second order 
weighted residual integral boundary layer equation, which in particular may be used to describe 
eddies in the troughs of the wavy bottom. We present numerical results which show that our model 
is qualitatively and quantitatively accurate in wide ranges of parameters, and we use the model to 
study some new phenomena, for instance the occurrence of a short wave instability (at least in a 
phenomenological sense) for laminar flows which does not exist over flat bottom. 

1 Introduction 

The gravity driven free surface flow of a viscous incompressible fluid down an inclined plate has various 
engineering applications, for instance in cooling and coating processes. For a flat bottom the problem, 
governed by the Navier-Stokes equations, is extensively studied experimentally, numerically and analyt- 
ically, see, e.g., [CD02] for a review. In particular it is well known that there exists a stationary solution 
with a parabolic velocity profile and a flat surface, the so called Nusselt solution, which is unstable to 
long waves if the Reynolds number exceeds a critical value Rent = 5/6 cot a, where a is the inclination 
angle [Ben57, Yih63]. However, the Navier-Stokes equations in combination with the free surface are 
hard to handle and one is often not interested in the flow field but only in, e.g., the film thickness F . 
Thus there has been much effort to derive model equations for the evolution of F . Because of the long 
wave character of the instability, length scales of free surface perturbations are large compared to the 
film thickness. Therefore a small parameter e can be introduced to scale downstream derivatives. By 
an asymptotic expansion approach a scalar evolution equation for F was derived in [Ben66] and later 
corrected in [Lin74] . However, this so called Benney equation has finite-time blow-up solutions even 
at moderate Reynolds numbers, see [PMP83]. Nevertheless, asymptotically it can be used to check the 
consistency of improved models, see [SRQM06]. 

Besides the reduction of the Navier-Stokes problem to a scalar equation for the film thickness F a 
hierarchy of less drastic reductions has been studied, starting with so called boundary layer equations, 
see again [CD02, Chapter 2], for instance. An important step was the derivation of an integral boundary 
layer equation (IBL) by Shkadov in [Shk67]. He used the averaging method of Karman-Pohlhausen 
which consists of taking a parabolic velocity profile like the stationary Nusselt solution as ansatz for 
the downstream velocity component U and integrating the streamwise momentum equation along the Z 
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coordinate perpendicular to the bottom. This yields a system of two evolution equations for F and the 
local flow rate Q = UdZ. 

Although the IBL reproduces various experimental observations like the existence of solitary waves 
it shows the following inaccuracies: 

1. The predicted critical Reynolds number differs from the exact value by a factor 5/6. 

2. The IBL is not consistent with the Benney equation. 

3. The assumed parabolic velocity profile does not fulfill the dynamic boundary condition at second 
order. 

The first problem follows from a linear stability analysis which yields Rcrit, ibl = cot a. For the second 
point one derives a scalar evolution equation for F from the IBL. This can be done by enslaving the flow 
rate Q to the film thickness F and expanding it in powers of e, which gives a scalar equation for dxF 
differing from the Benney equation already at order e, see [RQM98]. The third problem is due to the fact 
that the parabolic velocity profile has its maximum at the free surface which implies d:^ (F) = 0. 

Recently there has been much effort to overcome these problems. Along [RQM98,RQM00,SRQM06] 
a two-equation model for F and Q has been derived by a Galerkin method. Based again on a long wave 
expansion of the Navier-Stokes equations, the Nusselt solution and three more polynomials appearing in 
the derivation of the Benney equation served as ansatz and test functions. The resulting model consisted 
of four evolution equations for F, Q and two other quantities measuring the deviation from the paraboUc 
velocity profile. From this a simplified model, called weighted residual integral boundary layer equation 
(WRIBL) for F and Q was derived which is consistent with the Benney equation at order and pre- 
dicts the correct critical Reynolds number. However, this model does not reproduce well known solitary 
wave solutions if the Reynolds number exceeded a certain value only sUghtly larger than the instabiUty 
threshold. This deficiency can be cured by a Pade-like regularization method in [SRQM06]. Moreover, 
in numerical simulations the extension of the WRIBL to three-dimensional flows yields excellent agree- 
ment with recent experimental results from [PN03], see again [SRQM06]. See also [OGN08] for further 
detailed numerical studies of this model. 

The problem over wavy bottom is studied much less extensively. For experimental results we refer to 
[Poz88,VB02,WSA03,WLA05,VMHM05,AVB06,WBH+08]. On the theoretical side, [WA03,WLA05] 
give an expansion of Nusselt like stationary solutions in suitable small parameters and an analysis of 
their stability. In [Tri98, Tri04, Tri07a] the problem is studied numerically by simulations of both the 
full Navier-Stokes problem and model equations derived in a similar way as in [Shk67]. Moreover, a 
detailed numerical stability analysis based on the Navier-Stokes equations has been carried out [TriOVb]. 
In [DO07] a scalar Benney hke model has been derived and studied numerically, and in [HBAW09] 
an IBL over wavy bottom has been derived using Shkadov's method. Finally, using the method from 
[RQMOO, SRQM06] a first-order WRIBL has been derived and studied in great detail in [OH08]. 

Here we continue into a similar direction as [OH08] by deriving and analyzing numerically an alter- 
native WRIBL equation and a regularized version. However, in contrast to [OH08] our analysis is based 
on curviUnear coordinates from [WSA03] which allow to treat more general situations where for instance 
the free surface is not necessarily a graph over the (flat bottom) downstream coordinate. These curvi- 
linear coordinates are also more natural since they allow a clear distinction between flow components 
tangential and normal to the bottom. Moreover, our WRIBL is second order accurate which for instance 
allows the description of eddies in the troughs of the wavy bottom. Finally, our approach is somewhat 
simpler than the (more general) approach of [RQMOO, SRQM06] which consists of several polynomial 
ansatz and test functions in the Galerkin expansion. We find that by taking an accurate velocity profile 
U as single ansatz and test function in the Galerkin method the WRIBL can be obtained in one step. 

Thus, the outline is as follows: In Section 2 we present the governing equations in curvihnear coor- 
dinates. Since we focus on film flow over bottoms with long wave undulations we assume the bottom 
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Figure 1 : Sketcii of the geometry and the curvilinear coordinate system. 



steepness and the non-dimensional wave number to be of order e, < e <C 1, and expand all equations 
up to 0{e^). In Section 3 we derive an appropriate velocity profile serving as ansatz and test function 
used to derive our WRIBL by the Galerkin method in Section 4, and in Section 5 we check the con- 
sistency of the resulting WRIBL with the Benney equation over wavy bottoms. From the WRIBL we 
derive a regularized version called rWRIBL in Section 6 by removing second-order inertia terms which 
otherwise may lead to some unphysical behaviour. In Section 7 we finally give some numerical results. 
First, in §7.1, by comparison with available experimental and full Navier-Stokes numerical data we il- 
lustrate the accuracy of our rWRIBL over wide parameter regimes, including the occurrence of eddies. 
Second, in §7.2 we illustrate two new phenomena, namely that the bottom modulation may introduce a 
short wave instability (in a phenomenological sense) not present over flat bottom (except for rather ex- 
treme parameter ranges), and that and how the free surface may cease to be a graph over the (fiat bottom) 
downstream coordinate. A short summary is given in §7.3. 

2 Governing equations 

Figure 1 illustrates the inclined film problem with an undulated bottom h. The liquid is assumed incom- 
pressible and Newtonian, the Cartesian coordinate system , is incUned at an angle a with respect 
to the horizontal (a = 90° in Fig. 1), and the bottom profile h{x) is periodic with wavelength A and 
amplitude a. As we want to expand the governing equations in a small parameter e it is useful and nat- 
ural to introduce a curvilinear coordinate system for the following reasons. First, although the Nusselt 
solution is no longer a stationary solution if the bottom is undulated, for thin films and low Reynolds 
numbers the flow {u, w) is still mainly parallel to the bottom. To apply different scalings to u and w 
the coordinate system thus has to be orientated along the bottom profile such that the u component is 
tangential to the bottom, while using a fixed Cartesian coordinate system scaling involves a mixing of 
the Cartesian velocity components u, w. Second, for larger Reynolds numbers we may anticipate situa- 
tions as sketched in Fig. I where the free surface is not a graph over x and cannot easily be described in 
Cartesian coordinates. 

Thus, at every point of the bottom xe^ + h{x)ez we define a local coordinate system e^;, with e^; 
tangential and normal to the bottom. For an arbitrary point A within the liquid the arc length x of the 
bottom and the distance z along to the bottom are now taken as curvilinear coordinates. As we focus 
on fihn flow over weakly undulated bottoms this relation is always unique. Thus, 



in e^, coordinates, where 6 = 9{x) is the local inclination angle between and e^;. In order to 
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v(a;, z, t) = u{x, z, t)ex + 'w{x, z, t)ez 

f{x,t) 

p{x,z,t) 

Pair 

a 
P 

V 

g 



velocity field 

film thickness (perpendicular to the bottom) 

pressure 

pressure of the air above the liquid surface 
surface tension 
Uquid density 
kinematic viscosity 
gravity acceleration 



Table 1: Physical quantities. 



transform gradients we will also need the bottom curvature k which is defined by 

k(x) = ^ , ^. 

{i + {dib{x)r)-2 



(1) 



For further details concerning the transformation to curvilinear coordinates we refer to [WLA05]. 

To describe the free-surface flow we introduce the variables in Table 1. In contrast to Cartesian 
coordinates all quantities measured in curvilinear coordinates are written without a hat. The governing 
two-dimensional Navier-Stokes equations now read 



dtu + 



udxU + wdzU + 



1 + KZ 

1 1 



-KUW 



pi + KZ 
1 



1 + K,Z 

dxP + 9 sm{a — 9) + v 



1 



+ 



(1 + KZ)2 



{d^u - K^u + 2Kdxw) + 



(1 + KZ)^ 

1 



1 + KZ 



dxK{w — zdxu) 



udxW + wdzW — 



1 + KZ 



1 + KZ 



-KU 



= — dzP — g cos(a — 6) + ly 



(1 + KZ)^ 



dxK{u + zdxw) 



1 



+ J- — ^ — ^(^x^ ~ '^^^^ ~ 2Kdxu) + - — ^ — kOzW + d'^w 

( 1 ~1~ KZ] 1 ~\~ KZ 

-{dxU + Kw) + dzW = 0. 



1 + KZ 

At the bottom z = we have the no-slip and no-flux condition 

0. 



u\ 



z=0 



w 



z=0 



The dynamic boundary condition tangential and normal to the free surface z = f reads 



d-rW — KU 



0=((l + K/)2-a/)^ 

V 1 + 

(1 + Kf)d^xf - fdxfidxf - ((1 + Kf)^ + 2{dxff) K 
((l + ^/)2 + (a^/)2)3/2 

2pv 



a 



+ dzu]+ 4(1 + Kf)dxfdzw, 



+ (P - Pa 



l + {dxf/{l + Kf)y 



jdxffjdxU + Kw) dxf f dxW - KU 



(2) 



(3) 
(4) 

(5) 
(6) 



+ dzU 



(7) 
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while the kinematic boundary condition is 

d 



^^{f{x,t)-z) = ^ dtf + ^^^ud^f - w = 0. (8) 



In order to introduce dimensionless quantities we refer to the stationary solution over a flat incline. 

This so called Nuf 
thickness. We set 



f 2 

This so called Nusselt solution has the mean flow velocity {u) = ^^'3" — , where h is the constant film 



X='^x, Z=\z, F=\f, U=^u, 

A h h {y) 

^=;rT7T^' T=—^t, K=-^K, P=——p. 

Additional to a we can choose four non-dimensional parameters to write the governing equations dimen- 
sionless. To describe surface tension and viscosity effects we use 

Bi := — = ^ (inverse Bond number), 

A2 sin a pg\^ sin a 



{u)h 51/1'^ sin a 
"3^ 



R := = — —^5 — (Reynolds number). 



1 

Here lca= (^) ^ capillary length. The relation of Bi to the also frequently used Weber number 
W = — — is W = wBi. For the geometric quantities we introduce 

pp/i^ sin a o ^ J 

S := 277^ (dimensionless wave number), (, := 2Tr^ (bottom steepness). 
A A 

As we are interested in thin films over weakly undulated bottoms we suppose throughout that both 6 and 
C are of order e, where e is a small parameter, while R, Bi and a are assumed to be of order 1. The latter 
means that a is bounded away from zero such that cot(a) is bounded. However, a = 90° such that 
cot (a) = is allowed. 

All calculations will be exact of order e^, i.e. we keep all terms of order 1, (5, C, (5^, and 5C,. 
Throughout we will only display the C'(e^)-symbol if we want to emphasize that our calculations are 
only asymptotically correct. In all other cases we will skip it. In particular, skipping ©(e^) -terms, the 
dimensionless governing equations read 

SKdjU + 5Rd^U + 5Rd^W = -5RdxP + 3 ^^^f"~^^ + + dCKddJ + d^J, (9) 

sma 

S'^R&rW + S^RUdxW + S^RdzWW - SCRKU"^ = -RdzP - 3 ^"^f"~^^ + 6dlW, (10) 

sma 

dxU + dz{{l + 5CKZ)W) = 0, (11) 

;7(0) = W{0) = 0, (12) 

{l+2SCKF-S^{dxFf)dzU{F) + 5^dxW{F) - SCKU{F) + A6'^dxFdzW{F) = 0, (13) 

3Bi(aiF - ^K) = -R(P(F) - Pair) + 26dzW{F) + 0{6^), (14) 

drF+il- dCKF)dxFU{F) - W{F) = 0. (15) 



The dynamic boundary condition normal to the free surface (14), where we used the abbreviation ^ := j, 
is only given up to order e. As we are not interested in second-order terms of the pressure P this turns 
out to be sufficient. 
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3 A first-order velocity profile 

For given F we derive a solution [U, W, P) of the time dependent equations (9)-(14) which is exact 
to order e. By introducing the flow rate Q as independent quantity we also construct a velocity profile 
U which will serve as ansatz and test function in the Galerkin approach in Section 4. There, a first- 
order profile U = U{) + eU\\s sufficient since we can extract all necessary second-order terms from the 
boundary conditions. 

We assume that F is of order 1 while the velocity field {U, W) and the pressure P are enslaved by F 
and can be expanded in powers of e: 

U = Uo + eUi + 0{e^), W = Wq + eWi + 0{e^), P = Pq + ePi + 0{e^). (16) 

The geometric quantities K and 6 coming from the bottom profile can be expanded in powers of e, too. 
It turns out that the bottom curvature K does not contain terms of first order while the local inclination 
angle has a leading (, i.e. 

K = Ko + eK2 + oic''), e = Cdi + oie) 

with 9i{X) = dxB{X), see Appendix A. This yields 

— ^ = cot a (6*1 - zzC cot aOl + OiC^), — ^ = l-CcotaOi - zzC el + 0{C^). 

Since both 5 and C, are of order e, equations (9)-(14) read at 0{1) 

3 + d'pQ = 0, -RdzPo - 3 cot a = 0, d^o + dzWo = 0, 

C^o(O) = Wo(0) = 0, d^o{F) = Q, 3Bi{dj,F-^Ko) = -R{Po{F)-P^). 

The C'(l)-solution thus is 



Uo = -^Z^ + 3FZ, Wo = -^dxFZ\ Pq = | (cot a(F - Z) - Bi^iF + Bi^i^o) + Pair- 



(17) 



At 0{e) we get the equations 

SRdjUo + 6RdxUoUo + dRdd^oWo = -6RdxPo - 3Ccot a + ed'pi, 
- eRdzPi - 3Cei + 5dlWo = 0, 

d^i + dzWi = 0, Ui{0) = Wi{0) = 0, ddJiiF) = 0, -eRPi{F) + 2SdzWo{F) = 0, 
with solutions 



eUi = -5RdTF{Z'^-3F^Z) + 5RdxF -FZ^ - -F^Z 
2 Y 8 2 y 

+ 3(5 cot adxF - 5Bid\F + QBidxKQ + (cot aOi) {^Z'^-FZ ] , (18) 

eWi = - ^SROtxF - ^F^Z^^ + ^SROtFOxPFZ^ 

_ SRdj,F (^^FZ' - |P'Z') - SR{dxFf (^^Z^ - 3F^Z^ 

+ ^(5 cot adxF - SBidj^F + (BidxKo + Ccot a ei)dxFZ^ 

/I 3 

- {ScotadjcF - SBidjcF + (BidjcKo + CcotadxOi) ^Z^ - -FZ' 



2 2 



ePi = - ^Cei{Z -F)- ^ddxF{Z + F). 
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To get rid of the time derivatives of F we use the kinematic boundary condition (15) which leads at 0(1) 
to the identity 

OtF = -dxFUoiF) + Wo{F) + 0{e) = -SdxFF^ + 0{e). 
Thus Ui can be rewritten as 

- 3F''{ScotadxF-SBid^xF+C^idxKo+Ccotaei) {^-\ j • (19) 

If we assume temporarily that also the local flow rate Q = UdZ is enslaved by F we can easily 
state the £-expansion of Q = Qo + sQi + C){e^), namely 



Jo 



Qo = I UodZ = F^ (20) 

F 

ru dZ = . 
5 



/■^ 6 

eQi =£ / UidZ = -SRdxFF'' - F^((5cot adxF-SBid'^F+CBidxKo+Ccot a^i). (21) 

^0 



As mentioned in the introduction we cannot maintain the enslavement of Q to F since this would lead 
to a single evolution equation for F which fails to reproduce physics correctly. Therefore we treat Q as 
independent 0(l)-quantity and introduce a second representation 

U = U{F, Q) = Uo + eUi + ^(e^) (22) 

of the velocity profile which depends on both F and Q. For consistency, if we plug the enslaved version 
Q = Qo + eQi + C(e^) into (22) we must recover the expansion U = Uo + eUi + 0{e'^) calculated in 
(17), (19). This yields the following conditions for Uq: 

(i) / UodZ = Q as Q is of order 1 , (ii) Uq = Uq if Q = Qo + £)(£) is assumed. 
Jo 

As Q is independent of Z the first condition implies that Q occurs as a factor in Uq. From (20) we know 
that in the enslaved version of Q in zeroth order we have Q = F^. Thus 

3Q / 1 f ZV 



which is exactly the lubrication ansatz which is used in the method of Karman-Pohlhausen. Thus our 

new velocity profile will emerge as refinement of the parabolic profile. 
On the other hand, plugging Q = Qo + eQi into Uq yields 

- 3F\d cot adxF - dBid^xF + C^idxKo + C cot a ^i) + | j • (24) 

Thus, comparing (19) and (24), Uq contains terms which belong to Ui, and therefore Ui consists of less 
terms than Ui, namely 



-UI) +^(f) -if . (25) 
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To sum up, if Q is treated as independent 0(1) -quantity we obtain the first-order velocity profile 

Similarly, the second-order velocity profiles U2 and t/2 are derived in Appendix B. These are not needed 
for the derivation of the WRIBL but for the reconstruction of the flow field in Section 7. 



4 Galerkin method 

We start with the derivation of the evolution equation for F by integrating the continuity equation (11) 
along Z, i.e. 

/ djdJdZ + [(1 + 6CKZ)W]^ = 0. 

From Q = UdZ and the no-flux boundary condition we obtain dxQ-dxFU{F)+{l+5CKF)W{F) -- 
0, and eliminating W{F) by the kinematic boundary condition (15) and skipping all terms of order 
and higher finally gives 

9rF = -{l-SCKF)dxQ. (27) 

In order to derive an evolution equation for Q we first eliminate the pressure P from the streamwise 
momentum equation (9) before we apply a Galerkin method. By means of (10) P can be written as 



6RP{Z) = 6RP{F) -6R [ dzPdZ 

Jz 



= SRP{F) + 3j ^o^(" 6\dzW{F) - dzW{Z)). 

sm (y. 

To eliminate P{F) we use the dynamic boundary condition normal to the free surface (14) and the 
continuity equation (1 1) to obtain 

5RP(Z) = ^RPair + 5'^{dzW{F) + dzW{Z))-^B;{5d\F-(:K) + 3S ^^^["'~^\ f - Z) 

sm a 

= SRP^-6\d^iF) + d^{Z))-3Bii6dj,F-CK) + 3^ ^08(0-^) 

sma 

Plugging this into the streamwise momentum equation (9) we obtain 

6RdjU + 5RdxUU + 5RWd^ 

= 3?^^fc^ + 9lC/ + 2S'dj^ + 3SB,dl,F - 3CB,dxK - 3S'-^^^dxF 
sm a sm a 

- 35^^^^^^dxe{F-Z) + 5^^{dxU{F)) + SCKd^J. (28) 
sm a dX 

The next step is to perform a Galerkin method with the single test and ansatz function U from (26). 
Thus we plug U into (28), multiply the residual by U itself and integrate the result along Z. We want 
all calculations to be exact of order e^. This seems to be a problem since the first two terms on the 
right-hand side of (28) are of order 1 and we know U = Uq + eU\ + e^U2 only up to 0{e). However, 
the first term 3 is independent of Z, and by the definition of Q we get 

\'^^^^UdZ = 3'^^^^^Q. 
sin a sin a 
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The second term 9|C/ is slightly harder to manage. Integration by parts together with the no-slip condi- 
tion U{0) = yields 

nF f-F 

/ dyjjJdZ = ddf{F)U{F) - / {dzUfdZ, (29) 
Jo Jo 

and up to order the integral on the right-hand side reads 

= 3^ + ^/R'^{d^fQ' + Ge^l (1 - A) Q^,dZ. (30) 



At this point we need some information about the second-order term The velocity profile U 

emanates from the asymptoti 
UodZ = Q, which imph 
last integral in (30) satisfies 



emanates from the asymptotic solution U, and thus fulfills the boundary conditions (12), (13). Moreover, 
UodZ = Q, which implies U2dZ = 0. Therefore and due to the no-slip boundary condition the 



1 z 



^2 , d^2dZ 



F 1 rF 



which gives 



It remains to calculate the first term on the right-hand side of (29). From (13) we know that {F) = 
-5'^dxW{F) - 4.5'^dxFdzW[F) + 5C,KU{F) is of order where the velocity component W can be 
expressed by tj due to the continuity equation (11). Thus the C'(l)-terms of U are sufficient which means 
that we do not have to know U2 explicitly. This leads finally to 

The other terms in (28) are all at least of order e and we can calculate them rather easily by plugging in 

U = Uq + eUi. Testing (28) with U leads to 

SRdjQ = 5 sW«-g) ^_5Q_5^ cos(«-g) ^^^^_15^ sin(a-g) ^^^, 
2 sin a 2 2 sin a 16 sin a 

+ ^-B,{5d\F-CdxK)F-^^5R^dxQ + ^SR^dxF + ^J^^iQ 
+ %^CK§,+^^"§2(9xF?-66'p'xF-l6'^dxQdxF 

+ S'R'(^-^dxTQQF-^dTQdxQF-^{dxQfQ-^dlQQ^+^^dxQdxF^ 

where we made use of (27) to eliminate time derivatives of F. As there are still time derivatives of Q on 
the right-hand side this is not yet an explicit evolution equation for Q. However, from (20) we know that 
Q = F^ + 0{e), which leads to djQ = SF^drF + 0(e) = -3§dxQ + 0{e). Together with (27) this 
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gives the evolution system for (F, Q), namely 



dj,F = - {l-6CKF)dxQ, (31) 

SRdjQ = ^ ""("-^^ F - ^4 - 5^ cos(a-g) ^^^^ _ 15^ sin(«-g) 
2 sin a 2 2 sin a 16 sin a 

+ ^<5CK|+452^(axF)^-65^|aiF-^52^a^5xF-^5^R2(5^)^Q. (32) 

If we set the waviness = we obtain a system which is up to scaling the same as the non-regularized 
WRIBL in [SRQM06]. That means that in case of a flat bottom our one-step method is indeed equivalent 
to the Galerkin method with universal polynomials and subsequent simplification. Thus for C = our 
WRIBL is consistent with the Benney equation and predicts the correct critical Reynolds number Rcrit- 
In the next section we will check the consistency for > before we will regularize the equation in 
Section 6. 



5 Consistency 

The basic assumption throughout this paper is that F is of order 1 while U, W and P can be expressed 
in powers of e as stated in (16). In Section 3 this allowed us to solve the Navier-Stokes equations 
asymptotically, which was used in Section 4 to derive the evolution equation (27) for F depending on 
the flow rate Q. The natural approach to achieve a scalar equation is now to plug into (27) the expansion 

rF f-F f-F 

Q = Qo + sQi + e^Q2 + 0{e^)= U^dZ + e UidZ + U2dZ + 0{e^). 

Jo Jo Jo 

We call the resulting equation Benney equation for wavy bottoms. In (20) and (21) we have already 
calculated the zeroth and first order components Qq and Qi. Consistency now means the following: In 
the evolution equation (32) for Q we formally replace Q by an enslaved version Q^^^ with the expansion 



Q^'^^ = Qjf^ + £Qf ^ + e'Q^^^^ + 0(£^) . (33) 

It is remarkable that — is the only C'(l)-term in (32) which contains Q. Thus we obtain a set of 
linear algebraic equations for Qjf^, Qi^^, Q^^^ which can be solved easily. By plugging Q^^^ into (27) 
we obtain a second scalar evolution equation for F. We call our WRIBL consistent if this approach 
yields the Benney equation for wavy bottoms. 

To derive the Benney equation for wavy bottoms by a long wave expansion of the Navier-Stokes 
equations and the associated boundary conditions (9)-(14) we continue as in (17), (18). At 0{e'^) we 
obtain U2. As this is rather lengthy we refer to Appendix B and state here only the integrated version, 
namely 



£2Q2 = £2 / U2dZ 



I 

Jo 



^6'R'd],FF'^ + ^6'R\dxFfF^ + ^d^R{B,dj,F - cot adj,F)F' 

- —SQRiBidlKQ + cot adxOAF'^ + —5R(?,5B-AdlFf - 25 cot aidxFf 
35 5 

+h5Bid\FdxF-(:{B;dxKQ+ cot a ei)dxF)F^+^-^5'^RBi{dxFfd\FF^ 

5 

+ \5CKQF^-\5CdxeiF^+35'^d\FF^-5CeidxFF^+15'^{dxFfF^-\(^elF^. (34) 
8 8 2 
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Replacing Q in (31) by Qq + eQi + £^Q2 yields the Benney equation for wavy bottoms. 
Now we use (32) to derive a scalar model. Plugging (33) into (32) yields at 0{1): 

lF-l'^=O^Qr = F'- (35) 

At first order we get 

5 5 O^^^ 5 5 5 

SRdjQ^^ =--CcotaeiF--e^--dcotadxFF+-5Bid^xFF-^C^idxKoF 

- Y^R^dxQr + p^^-^dxF (36) 

By applying djQ^^^ = SF^drF = -3F'^dxQ^^+0{s) = -9dxFF^+0{£) this equation can be 
solved for Q^i^, which yields 

eQf L ^ (^SRdxFF^-C{BidxKo+ cot a 0i)-S cot a dxF+SBid^F^ F^. (37) 

Comparing these results with (20), (21) we already see that Q and Q^^^ match at zeroth and first order. 
In order to calculate Q^^^ we solve (32) at As this is somehow elaborate and does not give any 

new insight we state here only the result, i.e. = Q2 as expected. As both the long wave expansion 
and the WRIBL approach yield the same expansion of Q, the scalar evolution equations are in both cases 
the same. Therefore our WRIBL is consistent with the Benney equation also for C > 0. 



6 Regularization 

With the WRIBL (31), (32) we now have a second-order model for film flow over wavy bottoms which 
is consistent with the according Benney equation and reproduces in the limit of a flat incline the correct 
critical Reynolds number Rent- In order to achieve consistency the basic idea of the one-step Galerkin 
method was to use as test and ansatz function a velocity profile which is a solution of the expanded 
Navier-Stokes equations (9)-(14) also in the time dependent case. Therefore in (18) the first-order 
component Ui in particular contains the time derivative dxF which is substituted by the zeroth-order 
identity OtF = —SdxFF^. In contrast to setting dxF = in the velocity profile this procedure leads 
to the additional term —^d'^R'^{dxQ)^Q in the WRIBL (31), (32) which turned out to be necessary for 
consistency. 

However, over flat bottom it is known that a pure asymptotic expansion approach with the above 
substitution of dxF can lead to an unphysical behaviour if the Reynolds number exceeds a certain value 
Ro not far beyond Rcrit- In [PMP83] one-hump solitary wave solutions of a scalar Benney-like equation 
for flat inclines are considered. According to the bifurcation diagram [PMP83, Fig. 5] such homoclinic 
orbits are only found if the Reynolds number is close to the instability threshold, i.e. Rcrit < R < Ro- 
However, in [SAB94], where the two-dimensional Navier-Stokes equations were solved by a finite- 
element method, such a limit Rq was not obtained. Thus the asymptotic expansion equation used in 
[PMP83] appears to be valid only if R is not far beyond Rent, and shows non-physical behaviour if R 
exceeds a limiting value Rq. This deficiency appears to be closely related to finite-time blow-up solutions 
in the scalar Benney equation. 

For flat vertical walls it was shown in [SRQM06] using homoclinic continuation that such a limitation 
also occurs for the second-order WRIBL, i.e., the branch of homocUnic orbits again turns back if the 
Reynolds number becomes too large, see [SRQM06, Fig. 1]. However, if the inertia correction term, 
which corresponds to ~ i^S^R^ (OxQ)"^ Q in our notation, is neglected this non-physical loss of solitary 
waves ceases. At least for small C and otherwise similar parameters as in [SRQM06] we must expect 
similar problems with our model. 



11 



In [SRQM06] a Pade-like approximant technique is used to regularize the WRIBL in case of a flat 
incline, see also [Oos99] for the case of a scalar surface equation. The main idea is to remove the danger- 
ous second-order inertia terms by multiplying the residual equation for dtQ with a suitable regularization 
factor S. This procedure preserves the degree of consistency since the second-order inertia terms are still 
implicitly included. This becomes clear if one applies the zeroth-order identity Q = to S which 
yields the original non-regularized WRIBL. Homoclinic continuation now yields solitary wave solutions 
for the regularized model with no non-physical behaviour for R > Rent [SRQM06]. More precisely, 
for a wide regime of unstable Reynolds numbers solitary wave solutions are found, with amplitudes 
only slightly smaller than those obtained by numerics for the Navier-Stokes equations, in contrast to the 
regularization in [Oos99]. 

For the undulated bottom we again closely follow [SRQM06]. First, we spUt (32) into three parts, 
namely 

Resi:=5R{-dTQ-^^dxQ+^^dxF) and Res2 := -^{SR)\dj(QfQ (38) 
containing the inertia terms with leading SR and ((5R)^, respectively, and the rest 

2 sma 2 F'' 2 sma 16 sma 2 

The ^xQ-equation (32) now reads Reso + Resi -|- Res2 = 0, and using again Q = F^ v/e see that 
Res2 ~ {dxF)^F'^ is highly nonlinear. The aim is to get rid of the potentially dangerous term Res2 
without loosing the degree of consistency. Therefore, if we enslave again Q by F as in Section 5, 
no term up to 0{e'^) should be deleted or added. This is ensured, e.g., if we multiply the residual 
equation by a regularization factor S which can depend on F, Q and their derivatives. This yields 
5Reso +5'(Rcsi + Res2) = 0, and we are done if S fulfills 

^(Resi + Res2) = Resi +C'(e^). (39) 

This ansatz leads to the function 

S=il + ^Y- (40) 



Resi / 

Plugging the zeroth-order identity Q = F^ into (38) yields 

Resi = ^5RdxFF^ + ©(e^), Res2 = -^{5Rf{dxFfF'^ + 0{e^), 



and thus, using again Q = F 



3 



S ■.= (!- ^SRQdxF ) =S + 0{e'). (41) 



Then (39) leads to 5(Rcsi +Res2) = Resi +C'(e'^), and multiplying Reso + Resi +Res2 = by 6* 
finally yields the "regularized" equation 5Reso + Resi = 0(£^). In summary, the regularized version 
(rWRIBL) of the weighted residual integral boundary layer equation reads 

drF = - (l-5QKF)dxQ, (42) 

XT,/!^ 17xT,<3^^ , 9^_Q2 „ ^ /5sin(a-e)^ 5Q 5^cos(a-0)„ 
5Rd^ = - -SR-dxQ + -^5R-,dxF + [.—^F-^-^ - -^^—^dxFF 

lb sm a 2 2 lb r 

+^6'§-^{dxFr-65'p'xF-l6'^dxQdxF^ (l - ^SRQdxF^ . (43) 
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It is not easy to assess the value of this regularization. First, for flat bottom we numerically confirmed 
the loss of the one-hump solitary waves for the WREBL (31), (32) in a certain interval [Rq, Ri] of R > 
Rcrit and its regain for the rWRlBL (42), (43). However, here we use direct numerical simulations (see 
Section 7 for details), instead of homoclinic continuation in [SRQM06], which is not possible for C > 0, 
or in any case is much more involved since the solitary waves then do not decay to a constant state but 
to spatially periodic solutions. In these direct numerical simulations we find that the interval [Rq, Ri] is 
typically rather narrow, shrinks quickly with increasing C > and vanishes for ( greater some Co which 
depends on the other parameters. Also, the loss of solitary waves in [Rq, Ri] is not related to blow-up 
of solutions: instead, small amplitude irregular patterns appear in this interval. This might indicate a 
transition between two different branches of solitary waves for R < Rq and R > Ri, or some other more 
complicated structure in the background. 

To illustrate the effect of the regularization, Fig. 2 shows (in advance of §7) some differences between 
the rWRIBL and the WRIBL for a parameter set for which there is no interval [Rq, Ri] where the WRIBL 
does not have solitary wave solutions in direct numerical simulations. In general, these differences appear 
to be rather small, with the notable exception of the calculation of the critical Reynolds number Rent in 
Fig. 7 below, where the results for the rWRIBL are closer to available data. 

In general, in our simulations both the WRIBL and the rWRIBL did not show blow-up of solutions 
in parameter regimes of interest, but there appears to be one disadvantage of the rWRIBL: for some pa- 
rameters, as R becomes large the numerics for the rWRIBL fail more rapidly than those for the WREBL. 
In particular, for the parameters in Fig. 2 we can follow one-hump solitary waves for the WRIBL up to 
R « 90 where these split up into two humps, while for the rWRIBL we obtain numerical failures due to 
F — > pointwise for R not far beyond 12. However, this is strongly related to the method of simulation, 
i.e., to the fact that {F) = 1 is imposed, and should not be considered as blow-up of solutions of the 
rWRIBL: for instance we can follow one-hump solitary waves for the rWRIBL up to R = 21 if we 
double the domain length in Fig. 2. In summary, since we are more interested in the regime R not too far 
from Rcrit, where the rWRIBL gives results closer to available data than the WRIBL, below we focus on 
the rWRIBL for our numerical simulations. 




5 10 15 20 25 1 2 3 4 5 6 



(a) X [mm] (b) R 

Figure 2: Comparison of the WRIBL with the regularized version rWRIBL. a = 90°, 6 = 0.3, = 
0.05, Bi = 3.32 (comparable to [SRQM06, Fig. 1]); A = 5 mm, 5 bottom waves, R as indicated, and 
org and reg stand for the original WRIBL and the regularized version rWRIBL. (a) shows snapshots 
of the dimensionless film thickness F{x) with (F) = 1, and (b) the maximal amplitude of F extracted 
from one time period of well converged one-hump solitary waves. For these parameters, solitary waves 
of both the WRIBL and the rWRIBL are found for all R G (0.3, R2) with R2 f» 12, where (for the used 
discretization n = 400) the numerics fail for the rWRIBL due to F — > pointwise. Generically, the 
solitary waves for the rWRIBL have slightly smaller amplitude. 
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Figure 3: Local film thickness for two different inclination angles. For comparison with 
[WLA05, Fig. 3] it is measured not perpendicular to the bottom but to the main flow direction e^, see 
Fig. 1. Parameters: R = 0.0285, C = 0.31 and (a) a = 28°,S = 0.059, Bi = 2 x lO'^, (b) a = 18.05°, 
5 = 0.068,Bi = 3x 10-=^. 



7 Numerical simulations 

Though the rWRIBL (42), (43) is much simpler than the Navier-Stokes system (2)-(8), it is still a 
quasilinear parabolic system, with periodic coefficients. Therefore, a first step to explore some of its sta- 
tionary and non-stationary solutions are numerical simulations. For this we have set up a finite difference 
method with periodic boundary conditions in space for both, the rWRIBL and the WRIBL. To calculate 
stationary solutions {F, Q)s we use a Newton method starting at constant {F, Q) which corresponds to a 
Nusselt flow, which in contrast to the flat bottom case is not a stationary solution over wavy bottom. For 
the time dependent problem we may also use constant (F, Q) or perturbations of some {F, Q)s as initial 
data. We then use an implicit and adaptive time stepping. Depending on the flow characteristics, the 
spatial discretization was on the order of 50 (Fig. 7) to 400 (Fig. 12) points per bottom wave. Numerical 
convergence was checked by refining the discretization without perceivable differences in the solutions. 

7.1 Comparison with available data 

First we want to compare our results with available experimental and numerical data. Therefore we 
have to somewhat relax the assumption used in the derivation of the WRIBL that R and Bi are of order 
1 compared to (5 which are assumed to be small. However, similar relaxations often appear in the 
application of asymptotic expansions. In other words, one goal of the present section is to study how far 
the asymptotic expansion can take us. As said above, we focus on the rWRIBL since it gives shghtly 
better comparison with available data. 

We first simulate the stationary problem for fluid and geometry parameters taken from [WLA05], 
namely v = 1110mm^/s,/9 = 0.969 g/cm^, a = 20.4 mN/m. The bottom is a sine with wavelength 
A = 300 mm, amplitude a = 15 mm and trough and crest at x = 0, x = 150, respectively. Fig. 3 shows 
the resulting local film thickness which is the distance of the free surface to the bottom contour measured 
in e^-direction, see Fig. 1. As inclination angles we take (a) a = 28°, (b) a = 18.05°. Choosing the 
Reynolds number such that the maximum local film thickness is the same as in [WLA05, Fig. 3] we 
obtain stationary solutions {F, Q)s which for the film height are in perfect agreement with experimental 
data, see Fig. 4. 

In order to explore wider regimes of parameters and to get more detailed comparison also with full 
Navier-Stokes numerics we reconstruct the flow field using the second-order profile (47), (48) derived 
in Appendix B for sinusoidal bottoms. Following [Tri98], see also [Tri07a, Tri07b], we simulate the 
flow of liquid nitrogen over a vertical sinusoidal bottom with wavelength A = 1.57 mm and amplitude 
a = 0.0875mm. The fluid parameters axe u = 0.182mm^/s, /9 = 0.808 g/cm^ and a = 8.87 mN/m 
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which yield an inverse Bond number Bi = 17.92. As Reynolds numbers we choose R = 5 and R = 20. 
Again we achieve free surface profiles which are in good agreement with the Navier-Stokes numerics 
in [Tri98, Fig. 10], and also the flow fields are quaUtatively and semi-quantitatively reproduced correctly, 
see Fig. 5 and 6. Namely, there occurs a recirculation zone of correct size in the trough of the bottom 
contour if the Reynolds number is increased. 





Figure 5: Free surface and reconstructed flow field for stationary solutions of (42), (43) for (a) R = 5 
and (b) R = 20. The other parameters are 6 = 0.15 respectively 6 = 0.24, a = 90°, Bi = 17.92, C = 
0.35, A = 1.57 mm. 



Above we calculated stationary solutions (F, Q)s which, by analogy with the flat bottom case, must 
be expected to be unstable in the considered regime {a = 90°), see also [WA03, AVB06,Tri07a,Tri07b]. 
In the following we report on some numerical experiments to investigate the stability of stationary so- 
lutions and on some time dependent solutions in the unstable case. The standard approach to study the 
stability of {F, Q)s would be to calculate the spectrum of the linearization of (42), (43) around {F, Q)s, 
either numerically or analytically by expansion of first the stationary solution and then the eigenvalue 
problem in suitable small parameters. Eigenvalues of the linearization can then be calculated using Flo- 
quet theory. See [Tri07a] for a detailed parametric study of stability using this approach for an IBL, 
and [Tri07b] for the full Navier-Stokes problem. 

Here, since we are mainly interested in the shape of non-stationary bifurcated solutions in case of 
instabifity, to determine stabiUty of {F, Q)s we rather use a less systematic ad hoc approach. We numer- 
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Figure 6: Full Navier-Stokes numerics for the parameters used in Fig. 5. Reprint of [Tri98, Fig. 10], 
with permission from Elsevier. 



ically calculate (F, Q)s for various R, with fluid and geometry parameters fixed. Then, on a domain with 
eight bottom undulations, we apply a localized perturbation, let the system run, and determine stability 
by growth or decay of the perturbations. This yields a critical Reynolds number Rent in terms of the 
remaining parameters. 





A 


B 


C 


p[g/cw?] 


0.969 


0.969 


1.00 




24.1 


24.1 


1.00 


a[mN/m] 


20.0 


20.0 


70.0 


A [mm] 


108 


108 


10.0 


a[°] 


45 


10 


10 


Bi 


0.01 


0.04 


16.2 



Table 2: Parameters used to study stability of stationary solutions, with resulting inverse Bond numbers. 



Again we first focus on non-dimensional parameters from [WLA05], namely a = 45° and Bj = 0.01, 
using the dimensional parameter set A from Table 2, and calculate Rcrit as function of see Fig. 7. In 
agreement with [WLA05, Fig. 7], see also [AVB06, Tri07b], we find that the wavy bottom strongly 
increases Rent compared to the critical Reynolds number 5/6 cot a. over flat bottom. In particular, also 
the quantitative agreement with [WLA05, Fig. 7] is very good. Here the most notable difference between 
the WRIBL and the rWRffil occurs: Rcrit is somewhat larger for the WRIBL and hence the rWRIBL 
appears to be more accurate. 

Figure 8 shows time dependent solutions, with ( = 0.5 from Fig. 7, but for graphical reasons with 
only two bottom undulations. Over flat bottoms, for R > Rgnt the most prominent solutions are the 
(experimentally, numerically and analytically well known) traveUng pulse trains [CD02]. Also over 
wavy bottoms pulse hke surface waves develop, and the effect of the bottom waviness is a periodic 
modulation of the amplitude and speed of the pulses: (a) shows the decay of a locaUzed perturbation in 
the stable case, while (b) shows the emergence of a pulse in the unstable case. 

7.2 Some new predictions 

The numerics in §7.1 have shown that (42), (43) reproduces known phenomena qualitatively and quanti- 
tatively, in particular the appearance of eddies in troughs of the bottom for larger (, and the occurrence 
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1.2 
0.8 

0.2 0.4 0.6 

c 

Figure 7: Critical Reynolds number Rent as a function of the waviness ( for parameter set A from Table 
2. Along the critical values S varies from 6 = 0.035 (R = 5/6) to 5 = 0.048 (R = 2.2). [WLA05] 
denotes Rent from [WLA05], multiplied by 2/3 due to a different scaling. The critical Reynolds numbers 
were calculated with a tolerance of ±0.05. 




Figure 8: Numerical simulations in the sub- resp. supercritical case for parameter set A from Table 2 and 
C = 0.5 which gives Rent ~ 1-4, cf. Fig. 7; two bottom waves with periodic boundary conditions, (a), 
R = 1.1, f{x); for larger t the solution relaxes to a stationary solution, (b) R = 1.6, /(x); the solution 
is unstable and a traveling pulse evolves. 



of a long wave instability when the Reynolds number exceeds a critical value Rent as well as the increase 
of Rent with Next we consider a lower inclination angle for which we again investigate the stability 
of stationary solutions by the method specified above. Taking the same fluid parameters as in parameter 
set A but with a = 10° we get the critical values in Fig. 9 denoted by parameter set B. In contrast to 
Fig. 7 the critical Reynolds numbers are no longer increasing monotonously but reach a maximum at 
C ~ 0.17. For larger values of the bottom waviness Rerit decreases, and for > 0.23 it becomes less than 
the critical Reynolds number 5/6 cot a for flat bottom. 

Next we increase the inverse Bond number by choosing A = 10 mm and the fluid parameters of 
water, see parameter set C in Table 2 and the resulting critical values in Fig. 9. The dependence on ( 
turns out to be more pronounced than in Fig. 7. Figure 10 shows related time dependent solutions for 
some supercritical values. For small (, e.g. ( = 0.04 in (a), the instability is long wave (pulses), but 
for ( = 0.06 in (b) the perturbation evolves into a finite wavelength pattern. Thus, for a = 10° and ( 
larger than a critical value Co ^ 0.06 there appears a finite wave number instability, and the wave number 
increases as the bottom waviness becomes larger, see Fig. 10 (c)-(d). Since also the amplitudes of these 
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patterns are very small we conclude on a phenomenological basis that Fig. 10 (b)-(d) shows short wave 
instabilities, where, however, the following remarks apply. 

The linearization of (42), (43) around some (Q, F)s always has a Floquet exponent ^i(O) = from 
conservation of mass. In other words, /Ui(0) = since we have a family of stationary solutions {Q, F)s 
parameterized by the total mass M = /q F{1 + l5(KF)dX. IfK^ ii{K) is a parameterization of 
the Floquet exponents of the Unearization by wave number, then short wave instability in a strict sense 
means that unstable Floquet modes appear only in an interval ±K G {Ki,K2) with Ki > 0. However, a 
finite wave number instability may also be due to a side band (i.e. long wave) instability, that is, a branch 
IJ,{K) of unstable Floquet exponents with Re/n(i^) = C2K^ — c/^K* + 0{K^) with C2,C4 > 0, which 
up to order gives = ^J^^ as the most unstable wave number. To distinguish this from a short 

wave instability one should actually calculate the spectrum. However, we take (b)-(d) as strong hints 
for a short wave instability, since if (b)-(d) were due to side band instabilities we would expect larger 
amplitudes. In any case, to distinguish (a) from (b)-(d) we may call the latter short wave instabihties in 
a phenomenological sense. 

Finally, if = C(l) is the wave number of a pattern for the rWRIBL, then k = is the wave 

number in the dimensional Navier-Stokes system. Thus, if for instance h = 0(1) and a = 0(1) are 
fixed such that 6 = 2Trh/X = 0(e) and C = 27ra/A = 0{e) are small due to A = 0(e~^), then 
k = 0(e). However, even in this case, as already said at the start of §7.1, in applications we always have 
finite e. For instance, in Fig. 10 (b)-(d) we find fc = tt/IO, fc = 7r/8, k = 37r/20 [mm~^] as the basic 
wave numbers (the smallest possible wave number over a domain of length 80 mm being 7r/40 [mm~^]). 

Over flat bottom, short wave instabilities are only known for very small inclination angles, see [CD02, 
Section 2.3]. In particular, calculating the eigenvalues of the linearization of the rWRlBL (42), (43) 
around the Nusselt solution for the above parameters but C = by a Fourier ansatz we find no short wave 
instabiUty in case of a flat bottom. 
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Figure 9: Critical Reynolds number Rent as a function of the waviness C for parameter sets B (Bi = 0.04) 
and C (Bi = 16.2) from Table 2; for parameter set C the letters a-d indicate the values of ( used for the 
time dependent plots in Fig. 10. 

In Fig. 1 1 , for fixed R and varying ( we plot the minimal and maximal downstream velocities of some 
stationary solutions used in Fig. 9, which shows that these are continuations of the Nusselt solution. For 
C > Ci the minimal velocity Umin becomes negative which is an easy diagnostic for the existence of 
eddies. In particular, from (q < C,i we find that the short wave instability sets in before the appearance 
of eddies, which shows that the short wave instability is an effect of the wavy bottom on a Nusselt Uke 
laminar solution. Figure 1 1 (b) shows the stationary solution and reconstructed streamlines for the short 
wave unstable parameters C, = 0.4, R = 4.2. 

Finally, Fig. 12 illustrates a rather strongly unstable situation where due to a relatively large trav- 
eUng pulse the free surface is not a graph over x. This was one of the motivations to use curvilin- 
ear coordinates. Downstream the bottom maxima where the local incUnation angle is larger than 90° 
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Figure 10: Time dependent simulations for parameter set C, eight bottom waves, plots of the flow rate Q. 

(a) C = 0.04, R = 7.4 (long wave instability), (b) C = 0.06, R = 9.7 (short wave instability, four waves), 
(c) C = 0.08, R = 9.7 (short wave instability, five waves), (d) C, = 0.2, R = 6.1 (short wave instability, 
six waves). 




0.2 Ci 5 

(a) C (b) X [mm] 



Figure 11: Stationary solutions for parameter set C and R = 4.2. (a) Minimal and maximal downstream 
velocity ttmin, ^imax of Stationary solutions depending on For C, > C,i 0.38 eddies occur, (b) Free 
surface and reconstructed streamlines for = 0.4. 



(0 mm < X < 150 mm in Fig. 12) we find a bearing-out of the free surface as a pulse passes. This 
overhang is typically rather small since the pulse is small as it lost mass when it climbed "uphill" 
(150 mm < X < 300 mm in Fig. 12) to the maximum of the bottom. On the other hand, running "down- 
hill", the pulse grows and reaches maximum amplitude around x « 180 mm. This yields an overhang 
(to the left) of the free surface at the beginning of the "uphill" section. 

In the literature we did not find data or solutions comparable to Fig. 12, or to the short wave instability 
explained in Figures 9 to 11. Thus we think it will be interesting to study either experimentally or by full 
Navier-Stokes numerics the accuracy of these predictions. 
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Figure 12: a = 90°, A = 300mm, a = 20mm, hence C = 0.42, R = 10, 5 = 0.32; B; = 0.003 and 
initial data {F,Q) = (1,1). (a) Free surface over x, dashed line is the bottom contour, (b) Film thickness 
/ over X. 



7.3 Conclusions 

Using a Galerkin method with only one ansatz and test function we derived the WRIBL (31), (32) for 
film flow over wavy bottom, which in the limit of flat bottom equals the (one-dimensional version of the) 
WRIBL derived in [SRQM06]. In a second step we regularized the WRIBL to the rWRlBL (42), (43). 
Numerical simulations of the rWRlBL show very good agreement with available data from experiment 
and full scale Navier-Stokes numerics. Finally, our rWRlBL predicts two qualitatively new phenomena, 
namely a short wave instability of Nusselt like solutions (without eddies) at non-small inclination angles 
and at still rather small C, and solutions where the free surface is not a graph over the (Cartesian) down- 
stream coordinate. It remains to be seen whether these predictions can be verified experimentally or by 
full Navier-Stokes numerics. 

A Curvilinear coordinates 

In order to expand the non-dimensional curvature K and the local inclination angle in powers of we 
first scale the Cartesian coordinate x and the bottom profile b by 





. The relation between X and X is 
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thus X{X) = X - (dj^B{X)y dX + C>(C^), and therefore K{X) reads (cf. (1)) 

dpixix)) 



i+{d,b[^x{x))y 
1 



l + C^ (dj^BiXiX)))' 



X 



= -dlB{X) + -e\^3dlB{X){d^B{X)f + dlB{X) {dj^B{X)fdX j + 0{C') 
=: KoiX) + eK2iX) + O(C^). (44) 
For the local incUnation angle 9 we get 

9{X) = arctan (^dib (^-^X{X)j^ = arctaii(C5^5(X(X))) = CdxB{X) + O(C^) 

=:Cei{X) + 0{e). (45) 



B Second-order velocity profile 

Calculating the second-order component of the downstream velocity U = Uo + eUi + + C>{£^) by 
exactly the same approach as in Section 3 yields 

e^U2 = 6^R^dj,F (-ILfz^+—F^Z'--F^Z^+-F^Z^+-F^Z^--F^Z''+-F^Z 
^ \ 4480 560 20 40 8 10 7 



+5'K\dxFf 



.^Z^+ — FZ'~-F^Z''+-F^Z'+-F^Z^--F^Z^+—F'Z 
4480 560 80 10 8 5 35 



+S^R{Bidj,F-cotadj,F) (^^z'-^FZ'+^F^Z^-2F^Z^+^F'Z 

-5CR{B,d%Ko+cotadx0i)(^^Z'-^FZ'+lF^Z^-^F'Z'+^F'Z 
+6R {3SBi{djiFf-2Scota{dxFy+5SBid^xFdxF-C{BidxKa+cotaei)dxF) 
^FZ^-3F^Z^+6F'^Zj +6'^RBi{dxFfd'^xF r-Z^-l&FZ^+mF^Z 



+5CKo { ^Z^-^FZ^+3F^Z ) +5(0x01 



-\z'+\fz^-\f^z 



+5^dlF ( -Z^-^-FZ^+^F^Z^ +6CeidxF Qz2-3Fzj 



(^-^Z^+15FZ^ +C^el (^^Z^-'^FZ 



(46) 



If the flow rate Q is assumed to be enslaved by the film thickness F, then integration of (46) along 
Z G [0, F] gives the second-order component of Q = Qo + sQi + £^Q2 + 0{e^), see (34). In order 
to achieve an accurate velocity profile U2 depending on both F and Q we again treat Q as independent 
0(l)-quantity. This profile is not needed for the Galerkin method but only for reconstructing flow fields. 
Therefore we restrict our calculations to the practically relevant case of stationary flow over a sinusoidal 
bottom b{x) = a cos (^^x^ This impHes according to (44) and (45) 

K{X) = cosX + C»(C^), e{X) = -CsinX + C»(C^). 



21 



In case of the first-order profile (26) the correction of the parabolic profile turned out to be a self-similar 
polynomial with a coefficient depending on Q and dxQ but not on F or its spatial derivatives. The 
basic assumption now is that this is also true for the second-order correction IJ2. Therefore all spatial 
derivatives of F emanate from dxQ- As we consider here only stationary solutions the evolution equation 
for F (27) gives = 0. Thus in U2 we neglect all terms containing spatial derivatives of F. Taking 
again into account that treating Q as independent quantity mixes up e-orders in the expansion of U finally 
yields 



e^U2 = 5CR(cot a + B;) cos XQ 



+ 5C cos XQ 





3 
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- - 

4 




Thus the velocity profile used in Section 7 to reconstruct flow fields reads 




5Q cos XQ 




+ (5CR(cota + Bi)cosXg^ 40 ' ^ ' (l^ 



The according velocity component W is given by the continuity equation, i.e. 



(47) 



W 



l + 5CcosX Z 



f 

Jo 



dxUdZ. 



(48) 
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